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Abstract 

Space missions to small solar system bodies must deal with multiple perturbations acting on the 
spacecraft. These include strong perturbations from the gravity field and solar tide, but for small bodies 
the most important perturbations may arise from solar radiation pressure (SRP) acting on the spacecraft. 
Previous research has generally investigated the effect of the gravity field, solar tide, and SRP acting on 
a spacecraft trajectory about an asteroid in isolation and has not considered their joint effect. In this 
paper a more general theoretical discussion of the joint effects of these forces is given. 


1 Introduction 

Space missions to small solar system bodies such as asteroids and comets face significant dynamical chal- 
lenges. However these bodies are also of great interest to space scientists. Furthermore, small Near Earth 
Asteroids (NEA) are the most energetically accessible bodies in the solar system, and hence are a natural 
target for low cost space missions. The special challenges faced by spacecraft missions to these bodies are 
the relatively strong solar radiation pressure (SRP) forces that act on a spacecraft, the irregular mass distri- 
bution that can affect spacecraft in close proximity, and the gravitational effects of the sun, present through 
direct attraction and through the orbital motion of the asteroid in inertial space. In some mission scenarios, 
such as that encountered by the Japanese Hayabusa mission to asteroid Itokawa [1], these concerns lead 
to a mission design that does not have the spacecraft orbit the small body but instead maintains a quasi- 
stationary position in the asteroid-sun fixed frame enabled by frequent thruster firings. This introduces 
additional constraints, however, such as the need for frequent control thrusts, and also prevents precision 
estimation of the physical parameters of the asteroid in some cases. It also limits the lifetime of the orbiter, 
in that the spacecraft cannot stay indefinitely at the asteroid with no operator intervention. In the current 
paper we study the orbit mechanics of a spacecraft in orbit about a small body such as an asteroid. The 
purpose is to outline the limits that exist for stable orbital motion, so a clear determination can be made 
between when orbiting is feasible and when a spacecraft must instead hover about a small body. Specifically, 
we study a class of orbits that are stable against SRP perturbations over long time spans, and determine 
what the effect of the small body mass distribution is on these orbits. 

To analyze the stability of the proposed close proximity operations we apply a general analytical method- 
ology that rests on numerous studies of specific asteroids. We assume that the asteroid is modeled as a 
tri-axial body with 2nd degree and order gravity field coefficients, and incorporate the perturbative force of 
solar radiation pressure. It has been demonstrated previously that it is the 2nd degree and order gravity 
field terms that dominate the orbital stability for close proximity motion to uniformly rotating asteroids 
[2, 3, 4, 5], and that solar radiation pressure can destabilize spacecraft motion about a small body [6, 7, 8, 9]. 
The effect of both the gravity field and solar radiation pressure perturbations on the orbital evolution of 
a spacecraft can be characterized with analytical theories of motion. Using these theories we can develop 
a theory for characterizing whether close proximity operations about a small asteroid model is feasible or 
not. For characterizing stability against solar radiation pressure we apply the theory developed in [6, 7, 8] 
given the total mass of the asteroid, the asteroid orbit, and the mass to area ratio of the spacecraft. These 
provide us specific orbital limits and specific orbit designs to ensure a safe trajectory. For motion close to 
the body we apply theory developed in [2, 5] for the interaction of orbiters with non-spherical, rotating mass 
distributions. 
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Using the results reported here, questions of interaction between the solar radiation pressure and mass 
distribution perturbations can be taken into account. For a particular asteroid and spacecraft one can 
derive conditions under which there exists an interval of orbit radii that are robustly stable against both the 
destabilizing effects of solar radiation pressure (which becomes stronger for larger orbit radii) and interactions 
with the mass distribution (which becomes stronger for smaller orbit radii). These considerations lead to a 
nominal orbit design that places the orbit plane parallel to the asteroid terminator plane. This configuration 
is essentially a frozen orbit for motion about a small asteroid and yields constant orbit elements on average 
except for the longitude of ascending node, which becomes synchronized with the asteroid’s orbital motion 
about the sun. This synchronicity even persists for a highly eccentric asteroid orbit. When traveling through 
aphelion, however, the reduced precession rate of the spacecraft orbit can be more easily be perturbed by 
the precession rate induced by the asteroid mass distribution. This can lead to fluctuation of the orbit 
eccentricity and can make the orbit unstable. We develop specific bounds on fluctuations in eccentricity and 
discuss conditions for when this can destablize an orbit. 


2 Environment Models 

To begin our discussion we provide a brief review of the small body environment and force perturbations. 
These are split between the asteroid specification and the relevant force perturbations. 


2.1 Small Body Model 


We assume that the small body is on an elliptic orbit about the sun, and that its motion is well characterized 
by simple Keplerian dynamics over time spans of interest. This is generally true unless a close passage by a 
planet occurs, a situation that happens only infrequently. Of specific interest is the varying position vector 
between the asteroid and the sun, d = dd. Both the magnitude d, and a direction d are functions of the 
asteroid true anomaly v. 


1 + E cos v 

d = cos vx + sin uy ( 2 ) 


where x is the unit vector pointing to the orbit perihelion, y is the unit vector in the heliocentric plane of 

motion and normal to x. The cross product of these two vectors defines the orbit normal, specified as z, 

about which the asteroid revolves. 

Another important specification for the asteroid is its mass distribution and rotation pole. Here we 
assume the asteroid is uniformly rotating about its maximum moment of inertia. For our detailed analytical 
discussions we only consider the effect of the C20 = — J2 oblateness gravity field coefficient, although the 
effect of the C22 ellipticity gravity field coefficient is discussed in general terms. These gravity field coefficients 
can be specified by the mass-normalized moments of inertia of the body and equal: 

U 20 = ~( 2 I Z -I X -Iy) ( 3 ) 

C22 = \(Iy-I x ) ( 4 ) 

where I x A I y < I z are the principle moments of inertia of the body. 

The asteroid rotation pole is assumed to be fixed in inertial space. For convenience we specify its 
orientation relative to the above described orbital frame by an obliquity angle ft measured between the orbit 
normal z and asteroid rotation pole p, and a right ascension angle a measured between the vector z x p and 
x. With these definitions the asteroid rotation pole can be explicitly specified as: 


p = sin ft sin ax — sin ft cos ay + cos ftz 


( 5 ) 


2.2 Force Models 

Next we summarize the relevant force models that act on a small body orbiter. These include solar radiation 
pressure (SRP), mass distribution effects, and solar gravitational effects. 
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Solar radiation pressure: We use a minimal model for the SRP acceleration, modeling the spacecraft 

as a flat plate oriented at the sun. It is possible to generalize this model, however the essential dynamics of 
motion under this force will not be changed significantly. We start with a force potential formulation of the 
model: 


R g = gdr 


( 6 ) 


where g is the acceleration, d is the direction in which it acts (away from the sun) and r is the position of 
the spacecraft relative to the asteroid center of mass. The acceleration force is computed as the partial of 
R g with respect to r and is g = gd and the acceleration magnitude is computed as: 


9 = 


d 2 


(7) 


where /? = (1 + p)G\/B , G\ ~ 1 x 10 8 kg-km 3 / (s 2 -m 2 ), B is the spacecraft mass to area ratio in kg/m 2 , p is 
the reflectance of the spacecraft, /3 = (1 + p)Gi/B, and d is the distance between the sun and small body in 
km. This model assumes that the spacecraft is modeled as a flat plate facing the sun, and that the absorbed 
solar radiation is radiated away uniformly. Note the dependence of g on the distance between the asteroid 
and sun, thus the SRP force varies in time between a maximum at perihelion and a minimum at aphelion. 


Mass distribution: The attraction of the small body is initially modeled as a point mass, but we will 

later consider the effect of a 2nd degree and order gravity field perturbation. Inclusion of the 2nd degree and 
order field is sufficient to capture the main effects of non-sphericity in the mass distribution. Higher order 
gravity coefficients can be important, but do not have such a dominant effect on the qualitative dynamics 
of motion for this system. For our situation where the rotation axis of the body is not aligned with the z 
axis of the main coordinate frame, it is useful to state this term in vector form. 

R 2 O) = [! — 3(r - p) 2 ] +^C 2 2 [(f • s) 2 - (f • q) 2 ] (8) 

where we assume that the unit vectors p, q, s are aligned with the body’s maximum, intermediate, and 
minimum axes of inertia. For our analytical work we will neglect the C 22 coefficient, although this has been 
studied in detail elsewhere [5]. 

Solar gravitation: Also necessary to incorporate for some analyses is the perturbation of the sun’s gravity 

on the motion of an orbiter. This can be modeled as a 3rd body perturbation, and its functional form can be 
simplified by performing an appropriate expansion which is essentially Hill’s approximation. In this paper 
we do not deal with the solar gravitation in any substantial way, and we refer the reader to [7] for an in-depth 
discussion. 


3 Equations of Motion 


Combining the above force models we can define the equations of motion for a spacecraft in orbit about an 
asteroid. In an inertially fixed frame centered at the asteroid, they can be stated as: 


r 

b I sh 
II 

(9) 

u{ r) 

- - + R{ r) 

(10) 

R( r) 

= i? 2 (r) + R g {r) + R s (r) 

(11) 


We can equivalently state the problem in terms of the Lagrange Planetary Equations: 


da 

2 dR 


dt 

na da 


de 

1 - e 2 dR 

Vl^dR 

dt 

na 2 e da 

na 2 e du 


(12) 

(13) 
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1 


(14) 


di cot i dR 

dt na 2 y/l — e 2 du na 2 y/l 

duj Vl — e 2 dR cot i 
dt na 2 e de na 2 \J\ - e 

dQ 1 91? 

dt na 2 -\/l — e 2 sin i di 

da 2 91? 1 — e 2 91? 

dt na da na 2 e de 


dR 


— e 2 sin i 912 
dR 


(15) 

(16) 
(17) 


where a is the semi-major axis, e is the eccentricity, i is the inclination, u> is the argument of periapsis, 12 is 
the longitude of the ascending node and er = — M a where M a is the mean anomaly at epoch. The additional 
parameter n = \J fi/a 3 is the mean motion. 

For each of these equations, there are situations where we wish to pose them in a frame rotating with 
the sun-line d. By definition, this frame rotates at the true anomaly rotation rate, which varies with the 
distance of the asteroid with the sun. The rotational velocity vector for this rotating frame is z>z and 
v = \J fx aun A(l — E 2 ) / d 2 . For example, in a rotating frame the above equations of motion become: 


v r + 2vi ■ r r + i)z • r + v 2 i • z • r 


dU 

9r 


(18) 


where the subscript r denotes a time derivative with respect to the rotating frame, we note that v = 
—2g, sun E sin^/d 3 and that 


z = xy - yx 


(19) 


using a dyadic notation. The Lagrange Planetary equations are simply modified by subtracting the true 
anomaly v from the longitude of the ascending node, 12 r = 12 — v, yielding the modified equation: 


dl2 r 1 dR 

dt na 2 \J\ — e 2 sin i di 


(20) 


all else being the same. 


4 Orbit Mechanics 

The equations of motion stated above can be solved exactly for a number of different special cases. Some of 
these cases are of interest and can help answer certain fundamental questions such as what is the maximum 
orbit size before it is stripped away from the central body and what the secular evolution of an orbit subject 
to SRP. In the following we present some of these exact solutions and their applications. Much of this 
work has been reported elsewhere, although we bring it together for the first time and add some additional 
interpretation of these issues. 

4.1 Escape conditions 

The simplest system to analyze is that of a non-rotating, constant acceleration acting on a spacecraft orbiting 
about a point mass. This problem is integrable if motion is constrained to a plane containing the acceleration 
vector, as it is a limiting case of the fixed two-center problem with one of the centers being moved to infinity 
[10]. This problem was also analyzed extensively by Dankowicz [11, 12, 13, 14] using advanced methods 
of dynamical systems theory. Dankowicz’s analysis considered both the unperturbed problem and various 
perturbations of it. From the initial analysis of Dankowicz [11] we can extract some very useful results 
for spacecraft orbit design, which are based on the ideal conservation of the angular momentum about the 
constant force direction (i.e. , the line of action of the solar radiation pressure). 

The relevant force potential in this case is U = fi/r + gd ■ r where d is assumed to be stationary with 
respect to inertial space. The resulting equations of motion are: 

f r + gd (21) 


4 



We can easily show that the total angular momentum projected onto the d direction is conserved, or hd = 
d • (r x r) is a constant. The form of the equations are simplified if we shift to a cylindrical frame with the 
axis of the cylinder along d and the radius, p , and polar angle, 9, measured perpendicular to this direction. 
If we assign the x axis to the direction of the acceleration and the y and z axis perpendicular to this we find 
the simplified set of equations: 


X = 

px 

(22) 

1 

to 

II 

_pp 

r 3 

(23) 

p9 + 2 p() = 

0 

(24) 


where r 2 = x 2 + p 2 . The conserved angular momentum is now hd = p 2 9 and can be immediately inferred 
from integrating Eqn. 24. Eliminating 6 through this parameter we find the simplified set of equations and 
related force potential: 


x 

P 


dU d 

dx 

dU d 

dp 


1 A 


U d = --^+gx 
r 2 p z 


It is important to note that the reduced set of equations still have a Jacobi integral: 


C = \{x 2 + P 2 ) - U d 


(25) 

(26) 

(27) 


(28) 


This is directly related to the energy of the system, as if we substitute for hd 

C = \ (i 2 + P 2 + P 2 0 2 ) ~ U 


p 2 9 we find the following: 


(29) 


where U = p/r + gx is the original force perturbation and the velocity term x 2 + p 2 + p 2 9 2 = V 2 where Vj 
is the magnitude of the total inertial velocity. 

These equations have many interesting properties, but the one we focus on is the existence of an equilib- 
rium point that corresponds to a circular orbit, offset from the center of the point mass along the direction 
x and perpendicular to this same direction. A modified form of this orbit will also play a special role in 
the more general case accounting for the motion of the asteroid about the sun. The unique aspect of this 
solution is that it can lose its stability at a certain value of energy, and that this agrees well observed limits 
for the escape of a spacecraft from an asteroid due to the SRP perturbation. 

The equilibrium point is simple to find, solving for dUd/dx = 0 and dU d /dp = 0: 

(30) 

(31) 

where these two equations are still coupled through the relation r 2 = x 2 + p 2 . We note that the crucial 
parameter is the orbit radius, r, and the ratio g/ p. Once these are specified the value of x and p are fixed. 
We note that by definition x < r and thus we find a fundamental limit on the orbit radius, r < \fpjg. 
Conversely, as p < r as well, we find a lower limit as h 2 d / p < r, which implies a limit h d < 1/ v fpg. 

To study the stability of this equilibrium point we form the linearized equations of motion and compute 
the characteristic equation. Evaluation of the system yields a simple condition for stability: 


0 eq 


4 

P eq 


— r 3 

p ^ 

h A r z 

p eq 


V eq — 


3^3 V 


(32) 
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We note that the upper bound on x eq is the value at which the Jacobi integral value takes on its maximum 
value for all of the relative equilibria, C max = —2^/JIg/\/3. In [11] the orbit dynamics of the system were 
studied when the relative equilibrium were unstable, and escape was found to be the common occurrence. 
From the above analysis we see that when C < C max that there will be two relative equilibrium orbits, one 
stable with x < x* and one unstable with x > x* . When C = C max these two equilibrium points coincide 
at the maximum value of C . This simple relationship allows us to derive a necessary condition for escape of 
a spacecraft from the asteroid in this simple model, or a sufficient condition for stability. Simply put, if the 
Jacobi energy of a spacecraft is less than C max and the spacecraft is in the interior region of the zero- velocity 
surface, then it cannot escape. If its Jacobi energy is equal to or greater than C max , then it is possible for 
it to escape. In practice we find escape to be the usual situation. Thus, given an initial value of hd, x and 
p we have a limit for the reduced speed of the system, (x 2 + p 2 )/ 2 < C max + Ud- We can make this limit 
simpler and more useful by taking advantage of the special structure of this Jacobi integral, and find a limit 
on the osculating semi-major axis for stability: 


a < 



(33) 


This serves as a useful design parameter in constraining the maximum orbit size for mission design purposes. 

Once the motion of the asteroid about the sun is modeled, by allowing the unit vector d to rotate as a 
function of time, the conservation of hd is destroyed and the dynamics of the system become much more 
complex. Additionally, the solar gravitational effect must be taken into account as it balances the centripetal 
accelerations that arise from the frame rotation. In [7] this problem is derived and studied in detail. We 
will only summarize those results here. There are some exact solutions of the system in this case that 
are stationary in the frame rotating and pulsating with the asteroid-sun distance, and which correspond 
to periodic orbits in inertial space that follow close to the 2-body motion of the asteroid about the sun. 
Analysis of these equations provides a necessary condition for a spacecraft to escape which, in the limit of 
strong SRP perturbation, approaches the limit given above for the non-rotation case. Thus, this provides 
further validation of the above limit and makes it clear that it is appropriate to use for varying distances 
between the asteroid and sun. 


4.2 Averaged Orbit Mechanics Solutions 

Stating the spacecraft motion as a perturbation problem allows us to introduce averaging to the system. 
This allows one to solve for the general secular motion of an orbiter under SRP, which is found to be quite 
complex yet integrable, and to identify a set of “frozen orbits” that are suitable for orbiting a spacecraft in. 
We also present the averaged effect of mass distribution to identify the effect of the asteroid’s shape on an 
orbit. 

4.2.1 SRP Perturbation Formulation and Averaging 

Given that the perturbation from solar radiation pressure is usually small, it is useful to cast the problem 
into a perturbation form. The perturbing potential associated with SRP is found in Eqn. 6. We introduce 
the concept of averaging as this allows us to evaluate the secular effect of the perturbation on our system. 
The averaged potential is defined as 


1 r 2 * 

i? g = — y RgdM 


where M is the mean anomaly and e is the eccentricity vector and has magnitude equal to the eccentricity 
and points towards the orbit periapsis. The average is taken over the unperturbed 2-body motion of the 
spacecraft about the asteroid, and we hold the direction d fixed initially. 

Stated in this form, the rates of change of the orbit elements can be computed by substituting the 
averaged potential R g into the Lagrange Planetary equations. The simplest observation to make is that 


(34) 

(35) 
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the semi-major axis is conserved on average, and hence the energy of the orbit is conserved on average. 
The Lagrange Planetary Equations for this case can be integrated in closed form, a fact originally proven 
by Mignard and Henon [15], and worked out in detail in [6]. An alternate statement of the perturbation 
equations is given by [16] and uses the angular momentum and eccentricity vectors as the nominal orbit 
elements. These are not an immediately obvious choice, as this system is overdetermined in that the 6 
components of these two vectors are functions of only 4 orbit elements with secular rates, e, i, u> and O. 
Still, this formulation has a significant advantage in being able to be solved in closed form, a fact originally 
realized by Richter and Keller [16]. 

In [16] the standard set of equations for the secular evolution of the eccentricity vector e and the nor- 
malized angular momentum vector h in the presence of a non-rotating SRP force is found to be: 

6 = fvf c| - h (36) 

31 = < 371 

The normalized angular momentum vector h = \/\ — e 2 h, and thus is normalized by the factor yjjid, which 
is constant on average. 

These equations can be modified to account for the rotation of d. Given a frame rotational velocity 
vector of vi we find: 


e r + Oz ■ e 
h r + z>z • h 




3g 

— J-de 


(38) 

(39) 


where the r subscript indicates time derivative with respect to a rotating frame and will be dropped from 
this point on. We can write the secular equations as a linear system: 


e 

h 



(40) 


While this is a linear differential equation, it is not time-invariant as both v and g vary in time. It is 
important to note that the ratio i> / g is time invariant, as both vary inversely with d 2 : 


.[H = 

2^/nsunA{l - E 2 ) [ 

V a 

3/3 V 


3 g 

This quantity is a constant for a given asteroid, spacecraft and spacecraft orbit, so we define [15]: 


( 3/3 

tan A = — ,, . 

2 V HHSunA( 1 - E 2 ) 


(41) 


(42) 


We note that as the SRP perturbation becomes strong, A — > 7t/2, and that as it becomes weak A — + 0. 

Despite the time invariance of the ratio, the multiplying factor of the matrix is still time varying. We can 
eliminate this by changing our independent parameter from time to the true anomaly of the asteroid about 
the sun, as e = = e'v. Then the factor in front of the matrix becomes our newly defined SRP strength 

parameter, leading to the time-invariant linear differential equations: 


h' 


— z tan Ad 


e 

tan Ad — z 


h 


(43) 


where we note that the direction d is now fixed in the rotating coordinate frame. We find it convenient to 
redefine the independent variable from true anomaly v to a new angle ip = v / cos A. 
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(44) 


The solution matrix can be reduced to elementary functions, yielding the explicit solution: 


h (VO . 
$W0 


‘MV’ - ip 0 ) f° 

1A o 

cos {ip)I 6 xe + (1 - cos(V’)) 
cos 2 Azz + sin 2 Add 
— sin A cos A f zd + dz 

r 

- sin (ip) 


— sin A cos A ( zd + dz ' 
cos 2 Azz + sin 2 Add 


^zd + dz^j 


— cos Az sin Ad 
sin Ad — cos Az 


(45) 


This combines the explicit solution detailed in [15, 6] in classical orbit elements and the solution in [16] that 
only provides a solution in the non-rotating case. We note that the solutions are periodic in ip, repeating 
every true anomaly 27rcosA. Thus, over one heliocentric orbit the solution will advance 1/cosA times. As 
the perturbation grows large, and A approaches 7t/2, the solution will repeat many times. Conversely, as 
the perturbation grows small the solution will repeat only once every heliocentric orbit. Also, we note that 
the state transition matrix is ortho-normal and defines a rotation matrix in 6-dimensional space. 

We can discuss the existence of frozen orbits for this system. To find them we search for solutions to the 
algebraic equations: 

— z • e + tan Ad • h = 0 (46) 

— z • h + tan Ad • e = 0 (47) 


We note that e and h are orthogonal and assume that both are non-zero in magnitude initially. A detailed 
study of this matrix and its null spaces shows that there are two classes of solutions. First, if we choose e 
parallel to d and h parallel to z the second equation is identically solved and the first equation reduces to a 
vector equation parallel to the unit vector t. It turns out that the direction in which these vectors point is 
important, thus we introduce the test solutions e = e(d • e)d and h = h( z • h)z. Resolving the first equation 
along the direction t then yields: 

e(d • e) + tan A/i(z • h) = 0 (48) 

We note again that h = 1 — e 2 , which simplifies the expression to 

= -(d-e)(z-h)tanA (49) 

Thus there are two conditions for a frozen orbit to exist in this configuration: 

— (d • e)(z • h) = 1 (50) 

e = sin A (51) 

This class of frozen orbit were discussed in [6] and called Ecliptic frozen orbits. The orbit lies in the same 
plane as the asteroid’s heliocentric orbit. If the orbit normal is aligned with the heliocentric orbit, periapsis 
must point towards the sun, otherwise if the orbit normal is anti-parallel to the heliocentric orbit normal, 
periapsis must point away from the sun. As the perturbation strength grows the orbit approaches rectilinear, 
while if the perturbation strength vanishes the orbit approaches circular. Due to this, these orbits are not 
preferred for strongly perturbed situations, as the periapsis has a low altitude. Also, these orbits cross 
through the asteroid’s shadow. 

A second frozen orbit solution, called Solar Plane-of-Sky frozen orbits in [6], also exist and are more 
useful for highly perturbed systems. Now, choose e parallel to z and h parallel to d so the Equation 46 is 
identically solved and Equation 47 reduces to a vector equation again parallel to the unit vector t. Again 
the direction in which these vectors point is important, thus we introduce the test solutions e = e(z • e)z 
and h = h(d • h)d. Resolving Equation 47 along the direction t then yields: 


etanA(z • e) + h( d • h) 


0 


(52) 



which simplifies again to 


j x _ e 2 = -(z -e)(d • h)cotA (53) 

The two conditions for a frozen orbit to exist in this configuration are: 

— (z-e)(d-h) = 1 (54) 

e = cos A (55) 

The orbit lies in the plane perpendicular to the sun-line, commonly referred to as the terminator plane, a 
dawn-dusk orbit or a 3o’clock/9o’clock orbit. If the orbit normal points towards the sun, periapsis must 
point above the orbit plane along the positive z axis, otherwise if the orbit normal points away from the 
sun, periapsis must point below the orbit plane. As the perturbation strength grows the orbit becomes more 
circular, while if the perturbation strength vanishes the orbit approaches rectilinear. Due to this, these orbits 
are preferred for strongly perturbed situations. Also, these orbits avoid the asteroid’s shadow. Finally, these 
orbits are the natural continuation of the equilibrium orbits in the non-rotation case. 

4.3 Mass Distribution Perturbation and Averaging 

Now we will review the classical result for the averaged effect of a second degree and order gravity field. In 
order to relate this analysis to the SRP frame requires that the oblateness effect be specified in a general 
orbit frame, and not one chosen so that the inclination is measured from the symmetry axis. When a general 
frame for the oblateness perturbation not aligned with the symmetry axis is specified the inclination can 
also suffer a secular perturbation. To see this we state the oblateness perturbation function in a general 
vector expression. If the maximum moment of inertia axis of the mass distribution is stipulated as p and 
the potential is averaged, we find: 


R20 


L1C2 0 

2 a 3(l_ e 2 )3/2 


5-3 



(56) 


where h is the unit vector along the orbit normal. 

Now we wish to specify this force potential in our frame of choice, where we assume that the rotation 
pole of the asteroid, p, is specified by its obliquity angle f3 and its right ascension a and h = sin* sin fix — 
sin* cos fly + cos*z, so the general force potential in this frame becomes: 


R20 


MC20 

2a 3 (1 - e 2 ) 3 / 2 


5 — 3 (sin j3 sin * cos(a 


Cl) + cos (3 cos*) 2 


(57) 


Analysis of this potential shows that R 20 = constant, implying that h • p is constant which means that 
the orbit plane inclination relative to the rotation pole p is a fixed quantity. Also derivable from this form 
of the equations is that the precession rate about the orbit pole is a constant equal to 3nC 2 o/2p 2 h • p and 

that the argument of periapsis has a rate of advance equal to 3nC2o/4p 2 [l — 5(h • p) 


5 Stability of the Frozen Orbit Solutions 

5.1 Stability of the Frozen Orbit 

To study the stability of the frozen orbits defined above may seem to be a redundant exercise, given that 
the general solution of motion in these systems has been defined and is oscillatory. However, it is important 
to understand how these oscillations manifest themselves, especially when we consider the effect of joint 
perturbations between the SRP and mass distribution effects. For this stability analysis it is easier to work 
with the orbital elements themselves, as they do not have the non-linear constraints on magnitudes, etc., 
that the angular momentum and eccentricity vectors inherit. 
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We consider the two frozen orbit solutions in turn. In both cases we use the averaged force potential 
given in the form: 



e 


3 aeg 


[cos uj cos f l — sin uj sin f l cos *] x 
+ [cos oj sin + sin uj cos cos i] y + sin uj sin iz 


(58) 

(59) 


and the vector e is stated in our reference frame of choice where we align d = x . In the current discussion 
we do not shift to true anomaly as independent variable, as this will complicate the inclusion of the asteroid 
oblateness perturbation later. 

As a first step we restate the Lagrange equations for this perturbing force function, noting that a = 0: 


de 

3^_ 

dt 

2 na 

di 

3 3 

dt 

2 na 

du 

3 3 

dt 

2 na 

dfl 

3 3 

dt 

2 na 

da 

3g_l_ 

dt 

2 na 


: cos uj sin sin i 

= [(1 — e 2 ) cos uj cos S2 — sin uj sin f l cos i] 
■ sin uj sin O — is 


[cos uj cos f l — sin uj sin f l cos i] 


(60) 

(61) 

(62) 

(63) 

(64) 


where we note that this is stated in a rotating reference frame, captured by including the angular rate, is, of 
the asteroid about the sun. 

Due to its applicability, we only consider the stability of the terminator plane frozen orbit solution: 
i = 7t/ 2, sin il sin uj = —1, and e = cos A. Linearizing the Lagrange equations about this solution in terms of 
variables e, i, uj , and yields: 


dSe 

dt 


— sin AStt 

2 na 


d6i 

dt 


3 9 

2 na 


cot A Suj 


d6u 

dt 

d5fl 

dt 


3 9 1 

2 na sin A cos A 

33 1 r 

o — oe 

2 na sin A 


(65) 

( 66 ) 

(67) 

( 68 ) 


Thus, the system splits into two uncoupled harmonic oscillators. If we express these in terms of eccentricity 
and inclination with true anomaly as the independent parameter we find: 


5e" 

Si" 


1 


cos 2 A 
1 

cos 2 A 


Se 


Si 


(69) 

(70) 


with the oscillation period being 27 t cos A, which is to be expected given the general solution found for motion 
in this system. 


5.2 Perturbation from the Oblateness 

Now we want to consider the effect of central body oblateness on these frozen orbits. A reasonable assumption 
is that if the orbit semi-major axis is large, the effect of the asteroid shape may be small and negligible. 
When the asteroid is small, however, the maximum semi-major axis becomes small and must lie close to the 
body, raising the possibility for a destabilizing interaction between the SRP and oblateness perturbations. 
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We will evaluate this interaction by assuming that the spacecraft lies in a frozen orbit and then include the 
secular rates from oblateness as a constant perturbation in the Lagrange equations. The resulting system 
can be solved to find the maximum oscillation amplitude in the orbit element of interest, in this case the 
eccentricity and inclination. 

Our main interest is in the terminator frozen orbits, due to their favorable properties. Thus we will only 
study the effect of oblateness on this class of frozen orbit. To do so we must evaluate the rate of change of 
each of the orbit elements due to this perturbation at the frozen orbit conditions, stated previously. Doing 
so we find a = e — 0 and: 


di 

dt 

duj 

dt 

dll 

dt 

(±)n 


3 n(I a -I t ) , 2 


= (±) 


4 a 2 sin 4 A 

3 n(I a - I t ) 

4 a 2 sin 4 A 
3 n(I a — I t ) . 


4 a 2 sin 4 A 


sin (3 sin 2a 
[5 — 3 sin 2 / 3 sin 2 a] 
sin 2/3 sin a 


= sin 


(71) 

(72) 

(73) 

(74) 


We can insert these nominally constant terms on the right-hand side of Eqns. 65 - 68 to form a set of 
non-homogenous equations. Given the simple form of the linearized solutions, the general solution for these 
terms can be found assuming an initial orbit evaluated at the frozen orbit conditions. 

' Si 
Su> 

cos A {sin 2 (3 sin 2 a [sin(iz/ cos A) + 3 cos A(1 — cos(v/ cos A))] 

— 5 cos A(1 — cos{v/ cos A))} 

sin 2 (3 sin 2 a [(1 — cos(^/ cos A)) + 3 cos A sin (v/ cos A)] 

+5 cos A sin(z// cos A) 

and 

5e 

5Q 


. . 3 n(I a — It) sin 2/3 sin a 

(,-*-) 9’4a 

4 a 2 sm A v 

— cos 2 A( 1 — cos (v / cos A) ) 
cos A sin (v/ cos A) 


(76) 


3 n(I a - I t ) 1 

4 a 2 sin 4 A v 


(75) 


These results can be used to evaluate the fluctuation amplitude in these orbit elements from their nominal 
values, especially in eccentricity, in order to determine if the variation is acceptable. If the amplitude becomes 
too large the linearization assumption we have made will be violated, and a full non-linear evaluation should 
be made. We note that the amplitude of oscillation will change as the asteroid moves about the sun. The 
variation in eccentricity from its nominal value, cos A, is bounded by: 


|d'e| < 


3 n(I a — I t ) cos 2 A sin 2(3 sin a /2 
2 a 2 sin 4 A V HsunP 


(77) 


The deviation in e will decrease with an increasing semi-major axis, however there also exists an upper 
bound on a for the spacecraft to remain trapped. To determine the minimum deviation in eccentricity due 
to the asteroid’s oblateness, substitute the value a max = y/zjjJg/4: into the above to find: 


\Se\ 


mm 



{Ig - It) COS A 
H sin 3 A 


sin 2/3 sin a 


l 

d 2 


(78) 


where for this formula the assumption is that the orbit is adjusted to always be at the maximum distance. 
Thus, the second formula has the smallest perturbation at aphelion while Eqn. 77, for a fixed semi-major 
axis, has the maximum perturbation at aphelion. 
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5.3 Perturbation from Ellipticity 

For the ellipticity effect, we rely on the numerical analysis given in [5] where it is empirically determined 
that orbits outside of 1.5 resonance radii are not subject to destabiliziation due to the orbit ellipticity. The 
resonance radius is the orbit semi-major axis where the orbit period equals the rotation period and equals 
(/xT 2 /(47t 2 )) ' . Thus, for rapidly rotating bodies one can in principle orbit more closely, while for slowly 
rotating bodies one must in general maintain a greater distance. The 1.5 resonance radius limit we state does 
not make any assumption about the orientation of the rotation pole, and is only sharp if the orbit inclination 
is less than 45 degrees in general. For high inclination orbits, especially for retrograde orbits with inclinations 
greater than 135 degrees, it becomes possible to orbit much more closely to the body without suffering any 
effects from the ellipticity. In these situations, however, the oblateness becomes a major perturbation. Thus, 
the orbit limit to guard against ellipticity effects can be stated as: 


a > 


3 (T 2 fi\ 1/3 
2 ly UU ) 


(79) 


where T is the rotation period of the asteroid. In general, once the orbit pole of the asteroid is known it is 
possible to immediately map out when the terminator orbits will have to take special care relative to their 
interaction with the asteroid gravity field distribution. 


5.4 Destabiliziation Mechanism 

The mechanism for destabilization can include two different effects. First is the the oblateness alone, which 
may induce large oscillations in the frozen orbit elements which then excite the longer-term oscillations in 
the orbit elements which can lead to growth in the eccentricity. This form is not as serious, as it generally 
takes longer for the instability to become pronounced, and even then it follows a well defined solution. More 
difficult to deal with is the joint action of the oblateness and the ellipticity. Here, variations in the eccentricity 
can cause the orbit periapsis to drop, given that a is constant on average. If the amplitude becomes large 
enough, and the alignment of the orbit pole with the rotation pole is direct, resonant interactions between the 
spacecraft orbit and the asteroid rotation rate can cause changes in the orbit semi-major axis, eccentricity, 
and inclination. Variations in any of these can become reinforced, causing larger oscillations which send 
periapsis closer to the asteroid which further destabilizes the motion. When designing an orbit, this is the 
main mechanism to avoid, captured by the 1.5 resonance radius condition above, as this can destabilize an 
orbit rapidly. 


6 Example Computations 

In the following we provide some example computations for orbit mechanics of a spacecraft about two 
different asteroids, one relatively small and the other a bit larger. The spacecraft model we choose has a 
mass to area ratio of 33 kg/m 2 , corresponding to a projected area (including solar arrays) of 12 m 2 and a 
mass of ~ 400 kg. For simplicity we chose a reflectance p = 0. We modeled the motion of the spacecraft 
about two different asteroids modeled after NEA candidates for a mission. We call these Asteroid I and 
Asteroid II. 

Asteroid I has semi-major axes of 214 x 100 x 100 meters, an assumed density of 2 g/cm 3 , an assumed 
rotation period of 12 hours, and an assumed obliquity of 45°. It’s heliocentric orbit has a perihelion of 1.03 
AU and an aphelion of 2.7 AU. 

Asteroid II is larger and has semi-major axes of 635 x 317 x 317 meters, an assumed density of 2 g/cm 3 , an 
assumed rotation period of 19 hours, and an assumed obliquity of 45°. It’s heliocentric orbit has a perihelion 
of 1.1 AU and an aphelion of 1.45 AU. 
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Figure 1: Examples of an escaping orbit due to solar radiation pressure perturbations for Asteroid I. The 
left image is the trajectory as viewed from the sun, the right is as viewed in the terminator plane with the 
sun to the left. The spacecraft escapes when the asteroid approaches perihelion and a max is violated. 




Figure 2: Long-term stable orbits about Asteroid I (left) and Asteroid II (right). These orbits are viewed 
from the Sun to the asteroid and both are propagated for over a year, through the asteroid perihelion passage. 

7 Conclusions 

Spacecraft orbiting about small solar system bodies such as asteroids and comets must contend with signif- 
icant perturbations from solar radiation pressure, the body mass distribution, and solar gravitation. Orbit 
mechanics in the presence of each of these perturbations can be analyzed in detail, and in some special cases 
joint perturbation effects can be added. This paper states the known orbit mechanics results for motion in 
the presence of these perturbations and analyzes the effect of joint perturbations on the orbital motion of a 
spacecraft. 
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Figure 3: Two orbits about Asteroid II, one chosen in the terminator plane according to the nominal frozen 
orbit design and the other started 45 degrees out of the terminator plane in a circular orbit. The non-nominal 
design impacts on the asteroid in a little over 10 days. 



Figure 4: Example orbits about Asteroid II. Left, initial orbit radius of 2 km, at a distance where resonance 
effects from the mass distribution are important. Right, initial orbit radius of 3 km, outside the reach of the 
mass distribution effects. 
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